Phase Transition with Non-Thermodynamic States in Reversible Polymerization 
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We investigate a reversible polymerization process in which individual polymers aggregate and 
fragment at a rate proportional to their molecular weight. We find a nonequilibrium phase transition 
despite the fact that the dynamics are perfectly reversible. When the strength of the fragmentation 
process exceeds a critical threshold, the system reaches a thermodynamic steady state where the 
total number of polymers is proportional to the system size. The polymer length distribution has a 
sharp exponential tail in this case. When the strength of the fragmentation process falls below the 
critical threshold, the steady state becomes non-thermodynamic as the total number of polymers 
grows sub-linearly with the system size. Moreover, the length distribution has an algebraic tail and 
the characteristic exponent varies continuously with the fragmentation rate. 
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I. INTRODUCTION 

Equilibrium systems relax to a steady state described 
by the Gibbs distribution. In contrast, nonequilibrium 
systems are specified by the dynamics rather then by a 
Hamiltonian, and there is no general framework for de- 
scribing nonequilibrium steady states. Furthermore, un- 
like equilibrium phase transitions that are characterized 
by robust universality classes nonequilibrium phase 
transitions are highly sensitive to details of the underly- 
ing dynamics Q. 

In this paper, we investigate polymerization dynamics, 
and we report that competition between aggregation and 
fragmentation results in a remarkable non-equilibrium 
phase transition. Despite the fact that the dynamics are 
perfectly reversible, there is a nonequilibrium phase tran- 
sition from a thermodynamic state where the number of 
polymers is proportional to the system size into a non- 
thermodynamic state where the number of polymers is 
not proportional to the system size. 

Reversible polymerization is ubiquitous in polymer and 
atmospheric chemistry [3|, |J, [fj , and has analogies in net- 
works @ and computer science @, H, [1[ . Reversible poly- 
merization includes two competing processes: (i) The 
aggregation process [i] + [j] —* [i + j], merger of two 
polymer chains of lengths i and j into a larger polymer, 
occurs with the aggregation rate Kif, (ii) The fragmen- 
tation process [i + j] — ► [i] + [j] , breakage of a polymer 
into two smaller polymers, proceeds with rate This 
process is reversible because the aggregation process and 
the fragmentation process perfectly mirror each other. 

Reversible polymerization is described by the master 
equations [T(| 



dc k 
~~dt 



- K lj c l c J -c k ^2K t 



kj Cj 



i+j=k 



J>1 



2 * y 

i+j=k 



(1) 



where Cfc(t) is the density of polymer chains composed 
of k monomers at time t. The first two terms describe 
changes due to aggregation and the next two terms ac- 
count for changes due to fragmentation. The aggregation 
and fragmentation rates are (non-negative) symmetric 
matrices, = Kji and Fij — Fji. 

In the simplest case, the steady state distribution is 
found by equating the aggregation flux with the frag- 
mentation flux. 



K, 



F, 



(2) 



This detailed balance condition specifies an equilibrium 
state where the fluxes between any two microscopic states 
of the system balance. Such an equilibrium steady state 
exists for example when both th e ag gregation and the 
fragmentation rates are constant [ljj. Another equilib- 
rium state was found in a model of strings at very high- 
temperatures with the rates Kij = ij and Fij = i + j 
0. 

The detailed balance equation @ admits a solution 
only when the aggregation and fragmentation rates sat- 
isfy special relations, as shown in Appendix fX] In gen- 
eral, the steady state distribution is specified by the full 
master equations (|TJ) and moreover, the detailed balance 
relations |(5J) may very well be violated. For example, in a 
"chipping" process where only end-monomers can detach 
from the polymer, the matrix is sparse: F^ = when 
both i,j > 2. For constant aggregation rates, the chip- 
ping process exhibits a nonequilibrium phase transition. 
When the fragmentation rate falls below a certain thresh- 
old, a giant macroscopic polymer emerges (T3, ITil. [TH [lt| . 

We consider the aggregation and fragmentation rates 



K, 



X. 



(3) 



These rates, while intermediate between the linear chain 
model [llj and the string model [Hj], violate detailed 
balance (see Appendix [S]) . The product aggregation 
rate accounts for the natural situation in which any two 
monomers may form a chemical bond, thereby leading 
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to merger of their respective polymers. This polymeriza- 
tion process has been widely studied in polymer chem- 
istry [2ll ] and in the context of percolation 
[22I l23l 124 l25l] . The constant fragmentation rate reflects 
situations where all chemical bonds in the linear polymer 
are equally likely to break, thereby leading to breakage 
into two smaller polymers. This de-polymerization pro- 
cess has also been studied extensively [26[. Like the ag- 
gregation rate, the total fragmentation rate is linear in 
the molecular weight, Yl%+j=k ^ij = — 

Starting with N monomers, we study the nonequilib- 
rium steady states that emerge in the reversible poly- 
merization process (|3]). We find that the system gener- 
ally reaches a steady state, and that a nonequilibrium 
phase transition occurs at the critical fragmentation rate 
A c = 1. The average total number of polymers, -ZVtot, 
grows algebraically with the system size N, 



N 1 



with 



2A 
2 + A' 



(4) 



when fragmentation is weak, A < A c . The total number 
of polymers grows sub-linearly with the system size be- 
cause 7 < 1. Moreover, large polymers are likely as the 
polymer size distribution has a broad algebraic tail. The 
system develops this non-thermodynamic state through 
a gelation transition. We probe this gelation using mo- 
ments of the size distribution. 

In contrast, the system reaches an ordinary steady 
state when the fragmentation process is strong. The av- 
erage total number of polymers is proportional to the 
system size, iV to t = (1 — A _1 )./V, when A > A c . Large 
polymers become rare since the polymer size distribution 
has a sharp exponential tail. 

Interestingly even though the polymerization pro- 
cess is reversible because the underlying aggregation 
(M + \J] -* [i + j]) and fragmentation ([i + j] -> [i] + [j]) 
processes perfectly mirror each other and none of the 
transition rates ^ vanish, the breakdown of detailed 
balance leads to a remarkable phase transition involv- 
ing a non-thermodynamic phase where the number of 
polymers is not proportional to the system size and a 
thermodynamic phase where the number of polymers is 
proportional to the system size. 

The rest of this paper is organized as follows. The 
thermodynamic steady states that occur under strong 
fragmentation are examined in the next section, while 
the non-thermodynamic steady states that emerge when 
fragmentation is weak are analyzed in section IIIII The 
gelation transition is probed using the moments of the 
size distribution in section HVl Monte Carlo simulation 
results validating the theoretical predictions for the non- 
thermodynamic phase are detailed in section [V] We 
discuss the results and several open-ended questions in 
section I VII Appendices lATICl contain several technical 
derivations. 



II. THERMODYNAMIC PHASE 

Our focus is the steady state behavior and in particu- 
lar, the stationary polymer size density Cfc that satisfies 



- WW 



i+j=k 



fc Cfc = -A Cj + ^ (fc - l)c fc . (■")] 

j>k 



This steady state equation is obtained by substituting the 
aggregation and fragmentation rates ^ into the station- 
ary master equation ([TJ) . At the steady state, changes due 
to aggregation, represented on the left hand side, balance 
changes due to fragmentation, represented on the right 
hand side. Since both aggregation and fragmentation do 
not alter the total mass, the overall mass density £ fc fccfc 
is a conserved quantity, as follows from the rate equation 
([1]). We conveniently set the normalization £ fe fccfc = 1 
without loss of generality. 

The total polymer density M$ — £ fc Cfc is the most 
elementary probe for the state of the system. At the 
steady state, this quantity satisfies 



1 A , 

2 = 2^ 



M ). 



(6) 



an equation obtained by summing |(5]) and by using the 
identity £ fc J2j>k c 3 = £\ Cj £ fc<j 1 = £j(j - l)cj and 
the normalization condition £ fc fccfc = 1. The total den- 
sity is non-zero 



Mn = 1 - A" 



when the fragmentation rate is sufficiently strong, A > 1 . 
We focus on this strong fragmentation regime in the rest 
of this section. 

Let us assume that the system is large but finite with 
a total mass equal N, a state that can be achieved by 
starting with N monomers, for example. The expected 
total number of polymers, AT tot = N J2k °k = A^M is 
proportional to the system size N, and therefore, the 
system is in a thermodynamic state. 

The polymer size density can be calculated by utilizing 
the recurrent nature [27j of Equation For instance, 
the monomer and the dimer densities are 



Cl 



C2 



A- 1 
ATI' 
(A - 1)(3A • 



1) 



(A + l) 2 (3A + 4)' 



(8a) 
(8b) 



In general, the polymer size density is finite in the ther- 
modynamic phase, A > 1. Very large polymers are rare 
since the size distribution decays exponentially 



c fc ~ Ak~ 5 / 2 e- ak , 



(9) 



when k — > 00. This result is derived in Appendix [El 

When fragmentation is extremely strong, the system 
consists primarily of monomers: c\ = 1 + 0(A _1 ) and 
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c 2 = A -1 + 0(A~ 2 ) when A — > oo. The leading asymp- 
totic behavior can be obtained exactly in this strong frag- 
mentation limit, 



Cfc 



»fc-2 



k-1 



(10) 



for all k. This expression, obtained in Appendix [C] is 
compatible with the generic exponential tail @. 



The near-critical behavior 

The total density, the monomer density, and the dimer 
density all vanish near the transition point, Mo ~ (A — 1), 
ci ~ i(A— 1), and C2 — y(A— 1) as A — * 1. This behavior 
suggests the perturbative approach, 



Cfe = eb k 



(11) 



with the small parameter e = A — 1. The first two co- 
efficients are b\ — \ and b 2 = j. We substitute this 
form into the stationary equation §5§ and observe that 
the nonlinear aggregation term cx e 2 is negligible. Conse- 
quently, to leading order, the polymer size density obeys 
the linear recursion equations 



The size distribution is sharply suppressed according to 
© beyond this scale. Using the relation A ~ a -1 / 2 , 
derived in appendix [HI and a ~ fc^ 1 we deduce that 



Cfc 



-3/2^,-5/2 e ~constxe 3 k 



(17) 



for k 3> fc*. Indeed, this large size behavior matches the 
small size behavior ([15]) at the crossover scale 1)160. We 
conclude that the convolution term, that accounts for the 
creation of very large polymers from smaller polymers, is 
relevant only at very large scales, k 3> fe*. Otherwise, 
this term does not affect the density of small polymers. 

For completeness, we mention that the leading asymp- 
totic behavior of the moments, M n = J2kk n Ck, readily 
follows from the density fT5|) and the cutoff (|T5|) . 



M„ 



-3(r. 



n > 2/3 
n < 2/3. 



(18) 



Sufficiently large order moments diverge in the vicinity of 
the transition point, a consequence of the algebraic tail 
(|15p . The low order moments are finite, however. 



III. NON-THERMODYNAMIC PHASE 



kb k = ^ bj--(k-l)b k , 



(k + l)b k . 



/ j 3 
j=k+l 
oo 

-i = E ^ 

]=k+2 



2 kb k+ i. 



(12a) 
(12b) 



The second equation is obtained from the first by an in- 
dex shift. We subtract the two equations and obtain a 
recursion relation for the coefficients 



b k 



fc+i 
bk 



k + l 



(13) 



The coefficients can be conveniently expressed as a ratio 
of Gamma functions, b k oc T(k — 1/3) /T(k + 4/3), by 
using the identity T(x + 1)/T(x) = x. The polymer size 
density is therefore 



bk = - 



ir(i)r(fc- j) 
2r(|)r(fc. 



4)' 
3 I 



(14) 



where the proportionality constant is set by bi = h ■ 
Near criticality, the size density is algebraic, 

Cfc^efc" 5 / 3 , (15) 

over a substantial range, k <C fc*. This result follows from 
(fl4| and lim^^oo x a T(x)/T(x + a) = 1. Therefore, the 
likelihood of finding large polymers becomes substantial 
as the phase transition point is approached. The cutoff 
scale fc*, set by mass conservation, X)fc=i ^ c fe = 1) i s 
divergent 

k, - e~ 3 ■ (16) 



As the critical point is approached, the nonlinear con- 
volution term in ([5]) becomes irrelevant over the divergent 
scale (jTHJ) . By continuity, we deduce that the convolution 
term is negligible when A < 1. Consequently, the station- 
ary distribution obeys the linear equation 




(19) 



when A < 1. We introduce the normalized size density, 
Pk = Ck/J2k°k^ witn J2k>i Pk = 1- With this transfor- 
mation, the stationary equation (|19[) becomes 



k P k = A 1 



A 



(k - l)p k . 



(20) 



The monomer and dimer densities follows immediately, 



A 



Pi = 



1 + A' 



2A 



P2 



(1 + A)(4 + 3A) 



(21a) 
(21b) 



The normalized densities undergo a phase transition at 
A c = 1, as shown in figure 1. The fraction of monomers 
is not affected by the convolution term and (|21a|) holds 
for all A. However, the dimer density (|21b[) differs from 
the expression p 2 — (i+A^^+IjA) ^ or A < 1 implied by 
(|8bp and 1)7]). Similarly, the normalized size densities p k 
exhibit a phase transition for all k > 1. 
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FIG. 1: The normalized dimer density p2 versus the normal- 
ized monomer density pi . The phase transition at A c = 1 is re- 
flected by the discontinuity in the first derivative at p\ = 1/2. 



Generally, we recast Eq. (|20|) into the following recur- 
sion for the normalized densities 



Pk+l 
Pk 



U A_ 

fc 2+A 

, 2(1 + A) 
1 "r" 2+A 



(22) 



This recursion is obtained by repeating the steps 
leading to (p~5]> . Again, we express the nor- 
malized densities as a ratio of Gamma functions, 
Pk oc T(fc — 2^j)/r(fc + 2+A )■ ^he monomer density 
(|21a| sets the proportionality constant and hence, 



Pk 



a r(i+jg) r(fc-^ 
i + A r(i-^) r( fc+ ^) 



The size density has an algebraic tail, 



Pk 



k' 13 with /3 



2 + 3A 
2 + A ' 



(23) 



(24) 



as fc — » oo, thereby implying that large polymers are 
likely. The decay exponent l</3<5/3is not universal. 

The size density obeys c k oc p k , and the ^-dependent 
proportionality constant is obtained from the mass con- 
servation condition Y]u — i = 1 where the upper limit 
of integration is set by the system size. This sum is dom- 
inated by the density of large polymers. By performing 
the summation, we find that the polymer size density 
depends on the system size 



Ck 



(25) 



The total number of clusters N to t ~ Nj2 k Ck g rows sub- 
linearly with the system size N% Q t ~ N 7 with 7 = j3— 1 as 
announced in ^j. Therefore, the total polymer density, 
Mo ~ NP~ 2 , depends on the system size in contrast with 
the behavior when A > 1. In deriving the steady state 
equation (|19|) . we assumed that the convolution term is 
negligible. This assumption is consistent with the fact 



that the amplitude TV 3-2 in (|25|) vanishes as the system 
size diverges. 

Moreover, the expected total number of polymer of size 
k, C k = Nc k , is as follows, C k ~ N^k' p with 7 = 
This steady state is not thermodynamic! The number of 
polymers is much smaller then the system size, yet the 
number of polymers still diverges with the total mass: C k 
grows sub-linearly with the system size because 7 < 1. 

For irreversible polymerization, A = 0, all mass even- 
tually ends up in a single giant polymer, as reflected by 
the characteristic exponent 7 = 0. The total number of 
polymers still grows sub- linearly, N to t ~ -/V 2 / 3 , when the 
critical point is approached, A — ► 1. 

In the thermodynamic phase, all polymers are finite 
in size. Indeed, the exponential tail behavior ^ im- 
plies that the largest polymers are finite in scale. Near 
criticality, the scale of the largest polymers diverges ac- 
cording to (|16jl . In the non-thermodynamic phase, there 
are polymers of all possible scales because the power-law 
behavior (f2"5|) holds up to the system size oc N. Remark- 
ably, there are polymers that contain a finite fraction of 
all the mass in the system because according to (|2"5)l , the 
total number of macroscopic clusters, NJ2 k > CO nst xN c k 
is of the order one. 

The power-law distribution (f2~4"|) accounts for a compe- 
tition between two fluxes. There is a flux of mass from 
small scales to large scales that is generated by the aggre- 
gation process and a flux from large scale to small scales 
caused by fragmentation. The power-law behavior holds 
for all scales, indicating that these two fluxes balance 
at all intermediate scales. Similar competitions between 
the fluxes occur in fluid turbulence [281 ] . passive scalar 
advection [29[ , wave turbulence |30|, g ranular gases [3l[ , 
and driven aggregation systems |32l. l33l. l3~j] . However, re- 
versible polymerization differs from these driven system 
in that there is no external injection of mass to maintain 
the steady-state. 



IV. THE GELATION TRANSITION 

We now study the approach toward the steady state 
specified by the full master equation 

dc 1 00 A 

^ ^2 ijciCj - fcc fc + A^Cj - - (k - l)c fc . (26) 



df 



i+j=k 



j>k 



Initially, there are only monomers, c k (t = 0) = 5 k .i- 

First, we consider the total polymer density Mq that 
obeys the linear rate equation 



dt 2 2 1 ' 



(27) 



Subject to the initial condition Mq(0) = 1, the total den- 
sity is 



M = 1 - A" 1 +\~ 1 e- xt/2 . 



(28) 
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Hence, the steady state (JT)) is approached exponen- 
tially fast in the thermodynamic phase. Moreover, the 
monomer and the dimer densities also relax exponentially 
fast as in (f2"5| and in general, the polymer size density 
quickly approaches the steady state when A > 1. 

We focus on the kinetics in the more interesting non- 
thermodynamic phase where the total polymer density 
(|28p vanishes at time to = j In j^j- . This behavior is 
consistent with the vanishing total density, Mo ~ N 13 ^ 2 , 
implied by ([23)) . Of course, the negative expression l|28p 
is invalid beyond the time to- 

The moments M n = J2k ^™ c fc provide a direct probe of 
the kinetics. In the non-thermodynamic phase, the sys- 
tem nucleates large macroscopic gels and consequently, 
large moments diverge with the system size as follows 
from (|2"5|) . This, together with the vanishing overall den- 
sity Mo, indicates that the system undergoes a gelation 
transition at a finite time. At the gelation time t g , a gi- 
ant polymer or a gel emerges as is the case for irreversible 
polymerization (A — 0). From the master equation (|2"r?)) . 
the moments evolve according to 

dM n 1 ^— J / n A An— 1 

2 * — ' V m / 



2 n+ 1 



A A /n + T. 

—7 2^ I )SmAf„ + i- 



n + i * — ' \ to 



(29) 



where B m are the Bernoulli numbers [351 ] - For example, 
the second moment obeys dM n /dt = Mf — ^(M 3 — 1). We 
assume that large order moments diverge algebraically at 
the gelation time, M n ~ (i fl - for n > 1, and 

observe that the last term in the hierarchical equation 
(|2"9")l is negligible compared with the rest of the terms. 
We require that the time dependent term, the aggre- 
gation term, and the remaining fragmentation term are 
comparable and find that a = — b = 1. Hence, 

M n ~ (t g - t)-( n -V (30) 

for n > 1. Indeed, the moments diverge at a finite time. 
The exponent n— 1 characterizing this divergence is com- 
patible with the near critical behavior (TTHJ). The diver- 
gence (|3DJ} is different than the M n ~ (^ - ^)~( 2n - 3 ) 
behavior in irreversible fragmentation [23j, and there- 
fore, fragmentation quantitatively alters the nature of the 
gelation transition. 

We conclude that the solution to the master equation 
(|26p exhibits a finite time singularity. Even though we 
can not obtain the gelation time t g exactly, relevant prop- 
erties of the size density including the moments can still 
be obtained analytically. In particular, the form of the 
size density at the gelation time can be calculated by 
balancing the fluxes of mass due to aggregation and frag- 
mentation, following the scaling analysis in ref. [36j |. 

Consider M<"), the total mass density of polymers with 
size smaller then n, 



M ( - n \t) = Y J kc k {t). 



(31) 



According to the master equation (|26| . this quantity 
obeys 

dM^ 



dt 



E 



E 

i—lj-\-i—n — l 



-An(n + l)^ Cj . (32) 



j=n+l 

We now take the n — > oo limit. The aggregation loss term 
accounts for loss of finite size polymers to the infinitely 
large gel while the fragmentation term accounts for the 
balancing flux from the gel into small masses. We require 
that the two fluxes balance at the gelation point. We 
assume that at the gelation point, the size density decays 
algebraically 



Ck(t = tg) 



(33) 



for k 3> 1, as is the case for irreversible polymerization 
[24j . By dimensional counting, the aggregation flux term 
scales as n 5_2T while the fragmentation flux scales as 



,3-T 



The two fluxes balance when 5 — 2t 



t and 



as a result r = 2. Therefore, Cfe(t = t g ) ~ kr 2 , a behav- 
ior that is consistent with the aforementioned divergence 
of the moments (|3T)|) when the power-law behavior (|3"3"]) 
holds up to a cutoff scale that diverges near gelation, 
(tg-t)- 1 . 

The normalization condition J^ fe kc k = 1 imposes the 
exponent restriction r > 2 but the heuristic argument 
above yields precisely the marginal value r = 2. We 
therefore anticipate that there is a logarithmic correction 
with the following form ~ fc~ 2 (lnfc)~' 1 . We substitute 
this form into (|3"2"|) and then, the aggregation term is of 
the order n [lnTi] 1-2 ^ while the fragmentation term is of 
the order n [lnn]~ M . Therefore, fi = l, and [37j 

c fc - fc-^ln/c]- 1 (34) 

for k 3> 1. This decay is milder then the c k ~ k~ 5 / 2 
behavior found for irreversible polymerization [24|, and 
therefore, there are many more large clusters in reversible 
polymerization. 

The size of the largest gel at the gelation transition 
follows immediately from the extreme statistics criterion, 
NyifJLj Ck — 1. Remarkably, this size is nearly macro- 
scopic in the size of the system, 



iVtlniV] -1 . 



(35) 



fe=i 



This gel size is much larger compared with the k g ~ TV 2 / 3 
behavior in reversible polymerization. This increased 
scale enables the gel to withstand fragmentation. Wc 
also comment that the nearly macroscopic size scale (|35p 
provides the appropriate cutoff in (|3"Tj) . n ~ k g , and that 
at the gelation point, the upward mass aggregation flux 
and the downward mass fragmentation flux are nearly 
macroscopic, that is, they are proportional to the system 
size up to a logarithmic correction. 



The critical case 

For completeness, we discuss kinetics in the critical 
case A = 1 where the the cluster density (|2"5)l is purely 



6 



exponential, Mq = e - */ 2 . Similar decay characterizes the 
leading behavior of the size density. For example, the 
monomer density obeys dci/dt = — c\ + (M — ci) and 
consequently, 



Cl 



3 + 3 



2/ 



(36) 



relevant 



asymptotically, 
-t/a 



Only the first term is 

ci~|e~'/ 2 . In general, Ck ~ Uke~ t/2 , and by 
substituting this expression into the time dependent 
master equation (|26p . we observe that the nonlinear 
term is negligible. Consequently, the coefficients Uk 
satisfy the recursion equation 



2 



Uk 



j=k+i 



l)u fc . (37) 



From this recursion, the coefficients Uk satisfy Uk+i/uk = 
(fe — 2/3) /(k + 1). Therefore, the leading behavior of the 
size density is as follows 



Uk 



T(k 



3r(i)r(fe + i) 



(38) 



The tail of the size density matches the near critical be- 
havior (fTiaj). c k ~ e~ t/2 fc~ 5/3 . 



V. NUMERICAL SIMULATIONS 

We performed numerical simulations to validate the 
theoretical predictions. Below, we present results for the 
non-thermodynamic phase. 

The simulations were performed by starting with iV 
monomers and were carried by repeating the following 
Monte Carlo step. At each step, the total aggregation 
rate R a = s i s j an d the total fragmentation rate 

Rf = 2 X)i( s i~l) are calculated where Sj is the size of the 
ith polymer. Of course, both of these rates are propor- 
tional to the system size, R a <x Rf cx N. An aggregation 
event is executed with probability R a /{R a + Rf)i while a 
fragmentation event is executed with the complementary 
probability Rf/(R a + Rf). In each aggregation event, 
one polymer is chosen with probability proportional to 
its size and is merged with another polymer, also chosen 
with probability proportional to its size. In a fragmenta- 
tion event, a polymer, randomly chosen with probability 
proportional to the number of its bonds, is randomly split 
into two smaller polymers. Time is updated by the in- 
verse of the total rate t — > t + (R a + Rf) -1 after each 
Monte Carlo step. As a check, we successfully reproduced 
the total polymer density 0. 

In the non-thermodynamic phase, the system under- 
goes a gelation transition and then relaxes to the steady 
state. We run the simulations until the system relaxed to 
the steady state and then obtained the size distribution 
from a long series of measurements to reduce statistical 
fluctuations. The simulations results are for systems of 



size N — 10 5 . We present results for the normalized den- 
sities pk predicted in (|2"3")l . Overall, there is very good 
agreement between the theoretical predictions and the 
simulation results as shown in figure 2. The size distri- 
bution agrees with ([23]) and the tail of the distribution 
follows a power-law as in (JMJ). The simulation results 
agree with the theoretical results slightly better near the 
phase transition point (Figures 2a and 2b). Since the to- 
tal number of clusters grows sub- linearly with the system 
size, TVtot ~ N 1 with 7 < 1, extremely large systems are 
needed to reduce the magnitude of the statistical fluctu- 
ations. Such fluctuations are most pronounced at the tail 
region where the discrepancy between the theory and the 
simulation is a result of the limited system size. 

To quantify fluctuations in the total number of poly- 
mers, we also measured the variance er 2 = (N 2 ot ) — 
(Ntot) 2 - We find that the fluctuations follow a central- 
limit like behavior a 2 ~ (N tot ) (figure 3) and therefore, 



TV 7 



ith 



2A 
A + 2 



(39) 
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FIG. 2: The size distribution pk versus k for A = 0.9 (a) and 
A = 0.5 (b). The simulations results are for a system of size 
iV = 10 5 . 
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FIG. 3: Fluctuations in the total number of clusters. Shown 
are the average number of clusters, {N to t) and the variance, 



= (NL) - (iVtot) 2 for A = 0.5. 



VI. DISCUSSION 

In summary, we studied stationary and dynamical 
properties of reversible polymerization. We found an 
interesting phase transition involving a thermodynamic 
phase and a non-thermodynamic phase. When fragmen- 
tation is strong, the system is in a thermodynamic phase 
and the number of polymers is proportional to the sys- 
tem size. The system includes a large number of small 
polymers. When fragmentation is weak, the system is 
in a non-thermodynamic phase as the total number of 
polymers grows sub-linearly with the system size. In this 
phase, there is a small number of large polymers since 
the size distribution is power-law. Moreover, the poly- 
mer sizes are distributed at all scales. Macroscopic gels 
may exist as well. 

In the thermodynamic phase, the system quickly ap- 
proaches the steady state. The time-dependent behavior 
is much richer in the non-thermodynamic phase. The 
system exhibits a finite-time singularity: large moments 
of the size distribution diverge at the finite gelation time. 
At this time, the system nucleates macroscopic gels and 
the size distribution follows a universal algebraic decay. 
Past the gelation time, there is a second relaxation stage 
leading the system to a state where there are two bal- 



ancing fluxes of mass: aggregation transfers mass from 
small scales to large scales and fragmentation transfers 
mass from large scales to small scales. The stationary size 
distribution is algebraic but the characteristic exponent 
is not universal. 

Even though the aggregation process is described by 
non-linear terms, the analysis in the non-thermodynamic 
phase involves linear equations because the aggregation 
gain term is relevant only at the largest scale. Neverthe- 
less, formation of gels at the largest possible scale (in our 
case, the system size) is crucial in maintaining a station- 
ary state. Understanding the distribution of macroscopic 
gels is an open challenge, and our numerical simulations 
reveal an interesting anomaly with an enhancement of the 
population of macroscopic gels over the algebraic distri- 
bution (|25|) at the maximal scale. 

Missing from our calculations is a finite-size scaling 
analysis [38l [301. l40l l4l| near the phase transition point. 
Interestingly, the total number of clusters grows sub- 
linearly, Ntot ~ C(X) N 2 / 3 , just below the phase tran- 
sition point, but linearly just above the phase transition 
point, iVtot ~ (A — 1)N. Therefore, the amplitude C(A) 
may very well be divergent, and moreover, there must 
be an intermediate range of fragmentation rates centered 
around the critical point with a smooth crossover be- 
tween the two phases. The width of this transition region 
should vanish as the system size increases. 

Understanding fluctuations is another interesting di- 
rection for further research. We find Gaussian fluctua- 
tions in the thermodynamic phase and are able to com- 
pute the variance a 2 = jN but we are unable to obtain 
(|39[) in the non-thermodynamic phase. 

Finally, we checked that the same phase transition gen- 
erally holds as long as the aggregation rate is asymptoti- 
cally proportional to the molecular weights; for instance, 
when Kij = Aij + B(i + j) + C, a class of models that 
often arise in polymer chemistry 0, 0]. A more com- 
plete characterization of phase transitions in reversible 
polymerization is another avenue for future work. 
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APPENDIX A: DETAILED BALANCE 

In this appendix, we demonstrate that the detailed bal- 
ance equations © do not generally have a solution. The 
size densities Ck satisfy 



Knc\ = F n c 2 
K12C1C2 = Fi 2 c 3 
Ki3Cic 3 = F13C4 

K 2 2C\ = F22C4 



(Al) 



for k = 1,2,3,4. The dimer density and the trimer den- 
sity are uniquely expressed in terms of the monomer den- 
sity, 



K 



(-2 



11 „2 



cs = ^^4- (A2) 

-ril ^12 



However, there are two expressions for the 4-mer density 



c 4 



Kn K12 if 13 
Fn F12 Fj 



C4 



13 



(!rr) 2 lf c t- (A3) 



These two are identical only when the aggregation and 
the fragmentation rates satisfy the constraint 



if 12 if 13 _ if 11 if 22 

F12 F13 Fn F22 



(A4) 



Therefore, the detailed balance equations (|2|) have a solu- 
tion only for special aggregation and fragmentation rates. 
The constraint (|A4[) reflects the fact that there are mul- 
tiple paths between two states of the system. For ex- 
ample, a 4-mer may be formed by two dimers or by a 
trimer and a monomer. The detailed balance condition 
@ requires that the fluxes between any two states of 
the system balance along all possible paths. This con- 
dition leads to constraints of the type (|A4[) . The rates 
© violate the constraint (|A4[) as well as infinitely many 
other constraints. We conclude that in general, the de- 
tailed balance equations are overdetermined — there are 
infinitely many constraints like (|A4|) , and a solution does 
not necessarily exist. 



APPENDIX B: LARGE-SIZE ASYMPTOTICS 

We analyze the large-size asymptotic behavior of the 
polymer size-density in the thermodynamic phase by in- 
troducing the generating function 



G(z) = J2 



Ck e 



kz 



(Bl) 



The generating function obeys the differential equation 



(G 



IY2 



G' = A 



G - M e z _ X(G - G') 
1 -e z 2 



(B2) 
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where ' = 4z as follows from ([5]). Next, we shift the 
generating function by the total density 



G(z)=M + F(z), 



(B3) 



where Mo is given by ([7]) . With this transformation, the 
differential equation (|B2[) becomes 



(F 



r\2 



(2 + X)F' + 1 - A + XF 



e z + 1 



0. 



e* - 1 

By solving this quadratic equation for F' we find 

A 



F'(z) = 1 + - - y^{z) 



with the shorthand notation 



1 



$(z) = iA 2 + 2A-AF(z) 

4 e z — 1 



(B5) 



(B6) 



The asymptotic behavior of the size density follows 
from the singular behavior of the generating function. 
For instance, the asymptotic behavior 



c fc ~ Ak 



-a „ — ak 



(B7) 



implies that the generating function has the following 
expansion 

G(z) = G(a) + G'{a){z - a) + AT(1 - a)(o - z)^ 1 + ... 

(B8) 

when z — > a. Here, it is implicitly assumed that 
2 < a < 3. By differentiating this equation and by us- 
ing G' = F', we further obtain 



F'(z) = F'(a) + AT(2 - a) (a - z) 



a-1 



(B9) 



We now choose a to be the root of &(z), $(a) = 0, and 
as a result, equation (|B5j) becomes 

A 



F'{z) = l + --^¥{^{a~z) 1 l 2 + ... (BIO) 

when z — ► a. We obtain F'(a) = 1 + A/2 by matching 
the regular terms in (|B9|) and (|B 10|) . and 



= 5/2, ^r(-l/2) = - v ^¥ 7 (a) 



(Bll) 



by matching the singular terms. Therefore, the asymp- 
totic behavior is (|9|). 

The amplitude A can be expressed in terms of §'{a). 
Differentiation of equation (|B6[) yields 



*'(*) 



-A 



F'(z) 



1 



e'- — 1 

We next set z = a and use 



2F(z) 



F'(a) = l + -, 2F(a)^—tA = 4+-, 
2 ' , w e a - 1 2 



that follows from (|B6[) and $(aj = to obtain 



$'(a) 



(2 + A) 



1 



(8 + A) 



e 2a _ I 



(B12) 

(B13) 
(B14) 



By using equations (|B11|) and (|B14|) together with the 
identity T(-l/2) = -2T(l/2) = — s/4tt we obtain a rela- 
tion between the amplitude A and the constant a, 



/ A e° + l 
V 8?r e Q - 1 



2 + A - (8 + A) 



(B15) 



(B4) 

In particular, A ~ •^/3/(167ra), when A — > 1. 



APPENDIX C: EXTREMELY STRONG 
FRAGMENTATION 

The leading asymptotic behavior in the strong frag- 
mentation limit, A — > oo, can be obtained analytically. 
The steady state equation shows that c\ — 1 + 0(A _1 ) 
and C2 = A -1 + 0(A~ 2 ), and that in general, 



Cfc 



A— 1 



(CI) 



when A — > oo. 

To leading order, this form is consistent with the gov- 
erning equation ([5]) when the coefficients hk satisfy the 
recursion equation 



(k 



(C2) 



i+j=k 



The first two coefficients are /ii = 1 and = 1/2. 
We solve this recursion using the generating equation 
H(z) = Y^k kh k e kz . Next, we multiply JC2J) by fce fcz and 
sum over all k to find that the generating function satisfy 
the nonlinear differential equation 



H' — H = HH'. 



(C3) 



We now integrate this equation and find the implicit so- 
lution He~ H — e z . The explicit solution 



°° i.fc-i 



k=l 



k\ 



(C4) 



follows from the Lagrange inversion formula [42J . There- 
fore, the coefficients are hk = k k ~ 2 /k\ and the leading 
asymptotic behavior is (fTUjl . The large-size asymptotic 
behavior 



Cfc 



fc-i 



(C5) 



when fc 3> 1, obtained using the Stirling formula 
n\ r~i y/2Tmn n e~ n , is consistent with the generic asymp- 
totic behavior ([5]). 



